DEMO: HF calculations on small organic molecules by PySCF

DEMO: HF calculations on small organic molecules by PySCF#

We need to install a few things before we get started
- PySCF (for the quantum chemistry)
- py3DMol for visualizing the molecule
- plotly and kaleido for plotting

%pip install -q pyscf py3DMol plotly kaleido

Note: you may need to restart the kernel to use updated packages.

from pyscf import gto   # Used to define a molecule
from pyscf import scf   # Used to perform HF calculations
from pyscf import mp    # Used to perform Møller–Plesset PT calculations
from pyscf import cc    # Used to perfrom Coupled Cluster calculations
from pyscf import mcscf # Used to perform multireference calculations

Build and Visualize Molecules#

mol_xyz = """H      1.2194     -0.1652      2.1600
  C      0.6825     -0.0924      1.2087
  C     -0.7075     -0.0352      1.1973
  H     -1.2644     -0.0630      2.1393
  C     -1.3898      0.0572     -0.0114
  H     -2.4836      0.1021     -0.0204
  C     -0.6824      0.0925     -1.2088
  H     -1.2194      0.1652     -2.1599
  C      0.7075      0.0352     -1.1973
  H      1.2641      0.0628     -2.1395
  C      1.3899     -0.0572      0.0114
  H      2.4836     -0.1022      0.0205"""

mol = gto.M(atom=mol_xyz, basis="sto3g", verbose=3)

import py3Dmol

xyz_view = py3Dmol.view(width=400,height=400)
xyz_view.addModel(mol.tostring(format="xyz"),'xyz')
xyz_view.setStyle({'stick':{}, "sphere":{"radius":0.4}})
xyz_view.setBackgroundColor('0xeeeeee')
xyz_view.show()

#
# Use your mouse to interact with the molecule!
#

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Run calculations#

# Run HF and save the results
mf = scf.RHF(mol).run()

# In Jupyter notebooks you can hover over methods to see docstrings or you can print the docstrings out!
#print(mf.__doc__)

converged SCF energy = -227.890481261003

Visualize Results#

import plotly.express as px

# Plot the MO Occupations
fig = px.line(y=mf.mo_occ, markers=True, title="Molecular Orbital (MO) Occupations")
fig.update_layout(xaxis_title="Orbital Index (0-Based)", yaxis_title="MO Occupation")
fig.show()

# Plot the MO Energies (i.e. eigenvalues of the Fock matrix)
fig = px.line(y=mf.mo_energy, markers=True, title="Molecular Orbital (MO) Energies")
fig.update_layout(xaxis_title="Orbital Index (0-Based)", yaxis_title="MO Energies")
fig.show()

from pyscf.tools import cubegen

cubegen.orbital(mol, 'mo.cube', mf.mo_coeff[:,18]); # 42 electrons so 21 occupied orbitals

cube_view = py3Dmol.view(width=400,height=400)
cube_view.addVolumetricData(open("mo.cube").read(), "cube", {'isoval': -0.03, 'color': "red", 'opacity': 0.85})
cube_view.addVolumetricData(open("mo.cube").read(), "cube", {'isoval': 0.03, 'color': "blue", 'opacity': 0.85})
cube_view.addModel(mol.tostring(format="xyz"),'xyz')
cube_view.setStyle({'stick':{}, "sphere":{"radius":0.4}})
cube_view.setBackgroundColor('0xeeeeee')
cube_view.show()

#
# Use your mouse to interact with the molecule!
#

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Short Survey of more acurate methods#

from pyscf import gto, scf, dft, mp, cc, fci
import py3Dmol
import plotly.express as px

# Experimental geometry of gas-phase water
# Ref: https://cccbdb.nist.gov/expgeom2x.asp
mol_xyz = """O        0.0000   0.0000   0.1173
             H        0.0000   0.7572	 -0.4692
             H        0.0000  -0.7572	 -0.4692"""
mol = gto.M(
    atom=mol_xyz, 
    basis="6-31g", 
    verbose=4,
    charge=0,      # 0 by default
    spin=0,        # 0 by default, defined as (n_up - n_down)
    symmetry=True, # False by default
)

System: uname_result(system='Linux', node='fv-az1567-476', release='6.5.0-1025-azure', version='#26~22.04.1-Ubuntu SMP Thu Jul 11 22:33:04 UTC 2024', machine='x86_64', processor='x86_64')  Threads 4
Python 3.8.18 (default, Jul 16 2024, 19:04:00) 
[GCC 11.4.0]
numpy 1.24.4  scipy 1.10.1  h5py 3.11.0
Date: Fri Dec 13 17:45:35 2024
PySCF version 2.7.0
PySCF path  /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/pyscf

[CONFIG] conf_file None
[INPUT] verbose = 4
[INPUT] num. atoms = 3
[INPUT] num. electrons = 10
[INPUT] charge = 0
[INPUT] spin (= nelec alpha-beta = 2S) = 0
[INPUT] symmetry True subgroup None
[INPUT] Mole.unit = angstrom
[INPUT] Symbol           X                Y                Z      unit          X                Y                Z       unit  Magmom
[INPUT]  1 O      0.000000000000   0.000000000000   0.117300000000 AA    0.000000000000   0.000000000000   0.221664874411 Bohr   0.0
[INPUT]  2 H      0.000000000000   0.757200000000  -0.469200000000 AA    0.000000000000   1.430900621521  -0.886659497646 Bohr   0.0
[INPUT]  3 H      0.000000000000  -0.757200000000  -0.469200000000 AA    0.000000000000  -1.430900621521  -0.886659497646 Bohr   0.0

nuclear repulsion = 9.1895337629349

point group symmetry = C2v

symmetry origin: [0. 0. 0.]

symmetry axis x: [-1. -0. -0.]

symmetry axis y: [0. 1. 0.]

symmetry axis z: [0. 0. 1.]

num. orbitals of irrep A1 = 7

num. orbitals of irrep B1 = 2

num. orbitals of irrep B2 = 4

number of shells = 9

number of NR pGTOs = 30

number of NR cGTOs = 13

basis = 6-31g

ecp = {}

CPU time:        13.39

xyz_view = py3Dmol.view(width=400,height=400)
xyz_view.addModel(mol.tostring(format="xyz"),'xyz')
xyz_view.setStyle({'stick':{}, "sphere":{"radius":0.4}})
xyz_view.setBackgroundColor('0xeeeeee')
xyz_view.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Hartree-Fock #

Hatree-Fock (HF) is the starting point of the most of quantum chemistry
We variationally optimize the orbitals for a single Slater determinint
Working in the basis of atom-centered basis function we solve the Roothaan-Hall equations

See the PySCF user guide and examples for more info.

mymf = scf.RHF(mol).run()

Density Functional Theory #

In Density Functional Theory (DFT), the electron density of a reference noninteracting system is used to represent the density of the true interacting system.
The formulation resembles HF with a different effective Fock potential.
This effective potential depends on the density functional approximation which is chosen by the user.
PySCF gives users the access to a large number of functionals through the libxc and xcfun libraries.

See the PySCF user guide and examples for more info.

Møller–Plesset perturbation theory #

Perturbative corrections to the Hartree-Fock approximation.

See the PySCF user guide and examples for more info.

Coupled Cluster #

Perturbative method that improves on the Hartree-Fock approximation.
Coupled Cluster Singles and Doubles (CCSD) includes single and double excitation on top of the HF wave function.
Accuracy can be improved by including triples perturbatively (CCSD(T)).
Non-variational, but size extensive description of ground states. For excited states, see EOM-CCSD.

See the PySCF user guide and examples for more info.

mycc = cc.CCSD(mymf).run()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[12], line 1
----> 1 mycc = cc.CCSD(mymf).run()

NameError: name 'mymf' is not defined

e_ccsd_t = mycc.ccsd_t()

Full Configuration Interaction #

Full configuration interaction (FCI) is exact for a given choice of basis set.
Cost grows exponentially with the size of the system.
Also known as exact diagonalization.

See the PySCF examples for more details.

myfci = fci.FCI(mymf)
myfci.kernel();

Analysis#

Important data is saved to the PySCF method objects, making it easy to analyze and visualize!

# Collect data
methods = ["HF", "MP2", "CCSD", "CCSD(T)", "FCI"]
energies = [mymf.e_tot, mymp2.e_tot, mycc.e_tot, mycc.e_tot + e_ccsd_t, myfci.e_tot]

# Plotting
fig = px.line(x=methods, y=energies, title="Comparison of QC methods", markers=True)
fig.update_layout(xaxis_title="Method", yaxis_title="Energy (Ha)")
fig.update_traces(marker_size=12)
fig.show() # It's interactive!

DEMO: HF calculations on small organic molecules by PySCF

Contents

DEMO: HF calculations on small organic molecules by PySCF#

Build and Visualize Molecules#

Run calculations#

Visualize Results#

Short Survey of more acurate methods#

Hartree-Fock#

Density Functional Theory#

Møller–Plesset perturbation theory#

Coupled Cluster#

Full Configuration Interaction#

Analysis#

Hartree-Fock #

Density Functional Theory #

Møller–Plesset perturbation theory #

Coupled Cluster #

Full Configuration Interaction #